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Abstract 

Andreev bound states in clean, ballistic SNS and SNSNS junctions 
are calculated exactly and by using the Andreev approximation (AA). 
The AA appears to break down for junctions with transverse dimensions 
chosen such that the motion in the longitudinal direction is very slow. 
The doubly degenerate states typical for the traveling waves found in 
the AA are replaced by two standing waves in the exact treatment and 
the degeneracy is lifted. 

A multiple-scattering Green's function formalism is used, from which 
the states are found through the local density of states. The scattering 
by the interfaces in any layered system of ballistic normal metals and 
clean superconducting materials is taken into account exactly. The for- 
malism allows, in addition, for a self-consistent determination of the gap 
function. In the numerical calculations the pairing coupling constant for 
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aluminum is used. Various features of the proximity effect are shown. 

I. INTRODUCTION 

In 1991, Tanaka and Tsukada ||1| investigated the energy spectrum of the quasiparticle 
states in a superconducting superlattice based on a Kronig-Penney model potential. The 
calculations were done using the Andreev approximation |^ and the system considered was 
infinite in both the transverse and the longitudinal dimensions. In this paper we concentrate 
on three additional aspects: 2) we do the exact calculation and investigate the reliability of 
the Andreev approximation, ii) we do the calculation for a limited number of layers, which 
relates better to a possible experimental situation, and in) we investigate the transverse 
size dependence of the relevant properties such as the local density of states. The systems 
we consider cover the entire range from very narrow transverse dimensions to wider ones, 
thereby simulating 1, 2 and 3-dimensional systems. 

The authors are aware of two studies focusing on the transverse size dependence of 
Andreev bound states and the quasiparticle local density of states. Sipr and Gyorffy 
show the splitting of the degenerate bound states found in the Andreev approximation by 
studying the exact solution of the Bogoliubov-de Gennes equations. Their study is limited 
to a superconductor-normal metal-superconductor (SNS) junction and no self-consistency of 
the gap is considered. A second study, by Blaauboer et al. @], is a rather global one, but in 
that study the multiple-scattering Green's function formalism developed by Koperdraad 
was applied for the first time. The local density of states of normal metal-superconductor 
(NS) and SNS junctions were calculated, but the breakdown of the Andreev approximation 
was not noticed. 

We want to extend these studies by considering self- consistently obtained gap functions 
and by studying a system with more than two interfaces, looking for the influence of an addi- 
tional superconducting layer to the Andreev bound state spectrum. The multiple-scattering 
Green's function formalism to be used is an extension of the formalism introduced by 
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Tanaka and Tsukada [T[] in that no Andreev approximation is made. In addition it is set 
up in a quite different way, such that the multiple scattering of the (quasi)particles by the 
different interfaces is exhibited explicitly. 

The paper is organized as follows. First we give a concise account of the main features of 
the theory. In Sec. fT| the local density of states for the two systems considered is studied. 
In Sec. the gap function is calculated self-consistently for a bar-shaped superconductor 
and for the NS, SNS and SNSNS junctions. Concluding remarks are given in Sec. 0. 

Throughout the paper, Rydberg atomic units are used, in which the energy is in Rydberg, 
the distance is in Bohr (1 Bohr ^ 0.5 A), h = 1, and the electronic mass is |. 



II. THEORY 



Although in this paper a Green's function formalism is used, we first give the Bogoliubov- 
de Gennes equations 



E'^(r) = E 



( n(r) ^ 



f r 



(1) 



-V2-/i A(r) 
A*(r) V2 + /i 

used by Sipr and Gyorffy |Q to investigate the bound states of a SNS junction shown in the 
upper panels of FIG. ^. The main reason is that our formalism will be expressed in terms of 
the solutions of these equations. The inhomogeneity of the system is expressed by the space 
dependence of the gap A(r). In the superconducting parts the gap is a complex constant 
whereas in the normal metal it is zero. The spinor wavefunction describes quasiparticle 
excitations and the energy E is measured with respect to the Fermi energy /x. 

Application of the Bogoliubov-de Gennes equations to the bar-shaped superconductor 
shown in FIG. [I| yields 
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where (p is the phase of the complex constant A, and 
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ky = ^—;kz = — with Uy^Uz = 1,2,3,- ■ ■ (3) 

Ly Lz 



u''s = \ E + aJE^ -\A\2 (4) 



k^s = ik'F. + ^^E^-W (5) 
4. =/i-^?-/^' (6) 

The discretized nature of the transverse wavevectors is due to the vanishing of the wavefunc- 
tions on the transverse boundaries. The four standard solutions are labeled with the sign 
indices a and u that can both be equal to ±1. By this convention, the index a refers to the 
type of particle (electronlike for o" = +1 and holelike for a = —1) and the index u indicates 
the direction of propagation (z/ = +1 to the right and u = —1 to the left). The complete 
solution of Eq. (|l|) is a linear superposition of Eq. (H) for different possible combinations of 
(cr, z/). The above equations also hold for a normal metal by letting A = 0, in which case 
the subscript S is to be replaced by A^. 

The Green's function formalism is outlined extensively by Koperdraad et al. 0. It is an 
extension of the microscopic theory used by Tanaka and Tsukada, in that the electron-hole 
scattering properties are treated exactly, and it traces back to the microscopic description 
of superconductivity by Gor'kov and Ishii [0. For the sake of completeness and clarity 
we summarize its main features and add relevant new elaborations. The matrix Green's 
function satisfies the equation 



G(r,r',^u;„)=^(r-r')l (7) 



iun + + /i -A(r) 
-A*(r) iUn-V^-fi 

in which the differential operator is closely related to the operator in the Bogoliubov-de 
Gennes equations (|1|) apart from the replacement of E by iun where 

Un = nnksT with n a ± odd integer. (8) 

The quantity Un is called the Matsubara frequency. Possible inhomogeneities of the system 
are fully represented by the r dependence of the gap function. As far as the transverse direc- 



tions are concerned the general solution of Eq. (J^) can be expressed as a linear combination 
of solutions (0) over all allowed values of ky,k'y,kz, and k'^. Thus 
4 

G(r, r', iujn) = ^ G{x, x', ky, k'y, kz, k'^, iun) sin ky-y sin k'y-y' sin kzZ sin k'^z' . (9) 

J-'y-'-'Z h k' h h' 

The transverse standing waves are normalized to unity, which implies the orthogonality 
condition 

/■^^ - Ly 

sinkyysinkyydy = Yhy,ky^ (10) 

for the y coordinate and a similar condition for the z coordinate. The Fourier coefficient 
G{x, x', ky, k'y, kz, k'^, iuJn) can be obtained from Eq. using Eq. (plQl). It takes the form, 

4 r^v r^v r^^ r^^ 

G(x,x',ky,k' ,kz,k',iujn) = ^ ^ / sinfc^y / sin k' y' / sinfc^^; / sinfc'2;' 
^ LyLz Jo Jo ^ Jo Jo 

X G{r, r', iujn)dz'dzdy'dy. (11) 
Using Eq. (|TTp, Eq. can be written as 

G{x,x' ,ky,k'y,kz,k'^,iun) = 5(x - a;')4j„fc;4,,fc^l (12) 



^^n + £, + k% -A 



A* iuJn k%^ 

Before we proceed we want to make the terminology clear. G{v, r', iujn) is the actual Green's 
function. The Fourier coefficient G{x, x' , ky, k'y, kz, k'^, iun) is a Green's function in the sense 
that it is the solution of Eq. (|12D, but this equation demonstrates that it is diagonal in ky 



and kz- That facilitates both the notation, because the variables k'y and k'^ can be omitted, 
and the calculation of G{y,y' ,iujn) according to Eq. In calculating the local density of 
states and the self-consistent gap function, we will need the Green's function for diagonal 
spatial coordinates. As long as we keep x 7^ a;', we can already take y = y' and z = z' in 
Eq. (^. Since we are only interested in the variations over the longitudinal direction x, we 
can average over the transverse dimensions. By that Eq. (|^) simplifies to the series 

G{x,x',iUn) = G{x,x' ,ky,kz,iun)- (13) 



First we give the Green's function of a homogeneous bar-shaped superconductor 
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G'°(x,x',A;„A;„^a;„) =E^5^r(a^<)^r(a;>) = E4^r(a^>)^r(a;<), with 



(14) 



in which iVls = y {itOnY ~ I^P- The wavefunction ^pg'^ (x) is given by Eq. @ after omitting 
the transverse solutions, and and kg are given by Eqs. and (H) with the replacement 
E iujn- In addition the conjugate wavefunction is defined by 

Vs^i^) = ( M§e-*^/2 us'^e^'t'/^ ) e^'^^^'s^ (15) 



The Green's function ([T^ ) describes the propagation of excitations (holelike or electronlike) 
from the starting point x' to the final point x with the weighting factor d^. In this case of 
a superconducting bar, there is no scattering. 

Now we give the Green's function for a system with one interface. To account for the 
scattering at the interface, its appropriate form appears to be 

,ru-{x)t7^u'^ru'-W) (16) 

era-' 

where 

Glj{x, x', ky, kz, iuJn) = E dljVu'j{x)i)l,'f^{x') with (i = sgn(x - x'). (17) 

(T 

The label in Eq. (|14D has been replaced with the more flexible label uj indicating the 
position Xj of the interface. The subscript +j means x is in the right-hand side of the 
interface xj and —j means that it is in the left-hand side. The first term is nonzero only if 
the starting and final points, x' and x respectively, are at the same side of the interface. The 
second term takes into account the scattering of the particle at the interface. This scattering 
is aptly described by the parameter we call scattering t-matrix, tlj'jfj'. For a superconducting 
bar the t-matrix is zero, because there is pure propagation and no scattering. The scattering 
t-matrix can be obtained by imposing the continuity of the Green's function and its derivative 
at the interface. We then obtain 
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(18) 



In applying the above interface matching condition, it has come out to be convenient to use 
an extended definition 

^ o-„i<A/2 \ 
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(19) 



of the wavefunction ip'g' (x) which includes the derivative with respect to x. 

If there are more interfaces as in the systems depicted in FIG. multiple scattering 
occurs and the Green's function is given by 

Gujv'j'iX-i X , ky^ kzi iuJn) = G j^jj^,j,{x, X , ky, k^, iuJri){^vv'^jj' + ^~vu'^j+u,j'^ 

+ E (20) 

This is the generalization of Eq. ([T6|) . There are features in Eq. ( pO]) that do not appear in 
that equation. One of these features is the strange combination of Kronecker deltas in the 
first term. The first set of Kronecker deltas, S^u'^jj', insures that the first term is nonzero 
only if the starting and final positions are in the same layer. The other set, 6^u^>6j^uj>, serves 
the same purpose but at the same time it takes care of the redundancy in the indexing of 
the layers. The structure of the two sets of Kronecker deltas assures us that there is no 
overlapping of their functions. If one set gives the value unity the other set must be zero 
and vice versa. Another feature of Eq. (^) is the presence of the quantity T'^J^t^-t . Whereas 
the t-matrix describes the scattering of the particle at a particular interface, this quantity 
describes the multiple scattering of the particle at the interfaces along its path. We call this 
the multiple scattering T-matrix. The T-matrix is in fact a function of the t-matrix and is 
given by the multiple scattering equation 



rpaa'ufi' _ .acr'ufj.' fc X j_ X X \ i \ " laa" i/u" ja" rpa"(7' ,-u" fi' 



(21) 



To solve this equation it is necessary to first solve for the t-matrix at each interface using 
Eq. dg). 

Until now no approximations have been made in using the solutions of the Bogoliubov- 
de Gennes equations. In applying the formalism given above to the calculation of the density 
of states, iujn has to be replaced hj E + i6. After that one can make the frequently used 
Andreev approximation, which amounts to the replacement 

Je^ - lAP 

Kj - ^F. + (22) 

if k'^j occurs in the exponential and to k'^j ^ kp^ if k'^j occurs as a factor as shown explicitly 
in the third and fourth components of the wavefunction (p^Ql). This approximation is valid 
when -E, |A| << kp^. In the present paper we will investigate its limitation by looking at 
configurations in which E,A^ kp^. 

III. THE LOCAL DENSITY OF STATES 

Sec. provides the basic machinery to determine the matrix Green's function, which 
enables us to calculate the local density of states (LDOS) at a position x using the equation 

LDOS(x,E) = f^lim V lmGu{x,x;ky,k„E + i6) (23) 

in which Gu is the upper left matrix element of the multiple scattering Green's function (pO[). 
Although we only study Andreev bound states, which imply infinite peaks at the bound- 
state energies, the Green's function formalism makes it possible to broaden the peaks by 
using a small value of 5. The peaks acquire a finite height and must correspond to the 
bound-state energies. 

A. The SNS Junction 

We first apply the formalism discussed in Sec. ^ to a superconductor- normal metal- 
superconductor (SNS) junction. A schematic diagram of the system studied is given in 
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the upper and middle panels of FIG. 0. The dotted lines in the upper panel serve as an 
indication of the position of the inner NS and SN interfaces in the SNSNS system to be 
treated in Sec. [111B| . In this first application we will show more explicitly how the different 
labels are used. The interface index j has only two values 1 and 2. The T- matrix equation 
(pip can be recast into the form 



II II ^ — ^ / I II „ii „i ,,ii,,i 

-'--vjv'j' ~ ''-uju'j' ^jj' ' 2-^ ^-uju"'j'"'u"'j''' ~u" ,j+u" ,u''j'- 



(24) 



To implement this in matrix form, we must see to it that the elements of the T-matrix 
in both sides of the equation are arranged in the same manner. Thus, its elements in the 
right-hand side of the equation must be similarly displayed as in the left-hand side. This 



can be done by writing Eq. (|24D in the form 



rjiaa —vv rjiaa —vu 
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^-v2v'2 
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4.(7(7 —VV 



-v2v"2 ^v"2 



dZiioSu"- 



rpcra" — uv" rpaa" — uv" 



rpaa" — vv" rpaa" — vv" 
^-u2u"l -'-v2i/"2 



Each element of the above matrices is itself a 4 x 4 matrix with av (= ++, H — , — h. 



(25) 



as 



row and column indices. So actually, the matrices appearing in Eq. ( P^D have dimensions 
8x8. In this form the T-matrix can be obtained by appropriate matrix inversion. The 



Green's function is obtained using Eq. (|20D . The local density of states is determined by 
using Eq. (H). 

To investigate the validity of the Andreev approximation in our SNS junction, we will fo- 
cus on the choice of the transverse width Ly = = Lt of the junctions. The transverse wave 
components of the wavefunction (|]) are standing waves proportional to sm{kyy) sin{kzz) 
where ky and kz are given by Eq. (H). The different combinations of {ky, kz) or {ny,nz) are 
called modes whose allowed values are determined by 



TT 



^) inl + nl)>0. 



(26) 



When the transverse dimension is small, the second term in the right becomes large, as a 
result, only a few modes will be allowed. If this term exceeds the chemical potential /i, 
kp^ becomes imaginary, the wavefunction is damped, and consequently, such mode cannot 
exist. For larger transverse dimensions, the second term is smaller whereupon more modes 
are allowed. Most of our calculations will be done for small Lt so that only few modes will 
exist. We will tune Lt such that kp^^ is of the same order of magnitude as the gap energy 
A, in which regime the Andreev approximation (^2]) is not valid, and call such a Lt value a 
critical width. 

FIGs. D and ^ show the results for a configuration in which [uy, Uz) = (2, 2) is the highest 
allowed mode. The chemical potentials in the superconductor and in the normal metal, fis 
and fiN respectively, are assumed equal with magnitude 0.5. The longitudinal dimension 
L of the normal-metal part is 4000 Bohr and the gap A is treated as real with magnitude 
0.0001 Ry. The LDOS in the normal-metal part at x = 1000 Bohr is plotted against E/A. 
The peaks represent discrete energy states . We make the width of the curves, determined 
by the the parameter 6 in E + i6, wide enough so that the fundamental features can be seen. 
The numbers in parentheses denote the mode to which the energy belongs. In FIG. ^ the 
transverse width is determined by the condition that kp^ = A for the mode (2, 2), in which 
one finds that Lt = 12.5676 Bohr and in FIG. |^ the transverse width is Lt = 13 Bohr, which 
is slightly larger than the critical width but has the same allowed modes. In both figures the 
dashed curve is the local density of states calculated using Andreev approximation (AA) and 
the solid curves are the results from the exact calculation. In FIG. ^ the exact and the AA 
results coincide and just three states are found, one for each mode. For the critical width, 
the states for the first two modes are at almost the same position, but for the (2, 2) mode 
many states are found. We notice that the AA peaks have about the same height whereas 
the amplitudes of the exact peaks vary widely with the energy. The amplitude variation of 
the exact peaks is due to the specific position chosen for the LDOS. FIG. |^ depicts the exact 
LDOS at X = 1000 Bohr and x = 1500 Bohr. We observe that some of the peaks in the solid 
curve are suppressed while the corresponding peaks in the dashed curve are prominent and 
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vice versa. The suppressed peaks are pulled down by the small magnitude of the weighting 
wavefunctions at those values of x. The degenerate traveling waves corresponding to the 
AA are split in the exact treatment into odd and even (sine and cosine) functions having 
different heights at different positions 0. 

Finally, we want to comment on our choice of the transverse width up to now. In Sec. ^ 
it will come out that superconductivity is suppressed for transverse widths in the order of 
20 Bohr or less. This means that so far our choice of transverse widths may seem not too 
appropriate. We chose those transverse widths to illustrate with clarity the fundamental 
features of the bound states. In order to see the influence of a wider transverse dimension, 
we show in FIG. ^ the LDOS of a SNS system with = 100 Bohr (solid curve) and 
Lt = 99.8514 Bohr (dotted curve). The critical transverse width of 99.8514 Bohr is obtained 
from the condition kp^ = A for the mode (19, 12). To illustrate its main features clearly, 
we only show the results for the Andreev approximation. It appears that higher modes are 
allowed for a wider transverse dimension. The peak at about E = 0.95A comes from the 
many lower modes, each supporting just one bound state. Only the highest mode gives rise 
to the distinct set of peaks starting at about E = 0.074A. So effectively, the character of 
the results shown above is essentially unchanged. 

B. The SNSNS System 

A schematic representation of the system is shown in the lower panel of FIG. ^. The 
height of the middle superconductor is chosen to take the values hA where the range of h 
is < /i < 1. The LDOS is calculated using Eq. (^3]). The process of determining the 
t-matrix, the T-matrix, and the Green's function for the SNSNS junction is the same as 
in the SNS junction. The number of interfaces has increased by 2, and the dimension of 
the matrices has become 16 instead of 8. We again choose a square transverse cross section 
whose dimension is Lt = 12.5676 Bohr. This means that the set of modes is the same as in 
the SNS system. The lengths of the normal metal parts are each 1000 Bohr and that of the 
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middle superconductor is 2000 Bohr. 

FIG. shows the development of the bound states as the gap of the middle superconduc- 
tor is increased from 0.25A to A. All peaks belong to the highest (2, 2) mode, apart from a 
(1, 1) mode and a (1, 2) mode peaks just below i?/A = 1. The mode labels are shown only 
in FIG. ^ to avoid cluttering the figures. Both the exact solution and the states according 
to the Andreev approximation (AA) are shown. Again the lifting of degeneracy is observed 
in the exact results. The distribution of the peak heights shown corresponds to a calculation 
of the LDOS at x = —1500 Bohr. Results for other positions just give other distributions 
of peak heights. From now on, we concentrate on the position of the AA peaks to facilitate 
the comparison of the different figures. The SNSNS system with /i = is equivalent to the 
SNS system, so it is interesting to compare FIG. ^ with FIGs. 0a-d. One can see a general 
shift of the states in the latter figure to the right relative to the states in FIG. ^. For the 
system with h = 0.25 the h = state at E = 0.075A suffers a large displacement. This 
can be understood as follows. States below E = 0.25A see a longitudinal length Lj^ of 1000 
Bohr whereas the states above E = 0.25 A see a length of 4000 Bohr, L^r being the 
length of the N metal. According to a semiclassical result [§,0, bound-state energies scale 
as j^, so they are inversely proportional to the longitudinal length, which would imply a 
shift upwards to E = 0.3A. However, the separation of 2000 Bohr between the N part is 
smaller than the BCS coherence length ^ ^ 4500 Bohr, which suggests that the N parts 
of the system are not completely decoupled yet. This is the reason why the lowest state is 
found below E = 0.25 A, namely at E = 0.19 A. The positions of the other peaks are hardly 
changed. 

Looking at FIGs. |^-d one sees that the lowest state stabilizes at a position of about 
E = 0.26A. All other states are shifted more and more upwards for increasing h values. For 
h = 0.5, FIG. still a set of peaks is found above 0.5A, in FIG. |^c the set starts at about 
E = 0.75A. For h = 1, FIG. |^, only three states are seen just below E = A, apart from 
a state at E = 0.78A, lying three times as high as the lowest state, in agreement with the 
semiclassical picture. A test calculation for a SNS system with L^r = 1000 Bohr differs from 
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FIG. 1^ only as far as the position of the highest peak, the (1, 1) peak, is concerned. 

The idea of decouphng can also be illustrated by FIG. |^. For the sake of clarity we just 
give the results for the Andreev approximation. FIG. |^(a) shows the LDOS for a SNSNS 
system with L = 6000 Bohr. The lengths of the N parts are kept at 1000 Bohr but now 
the length of the middle superconductor is 4000 Bohr, slightly below the BCS coherence 
length of approximately 4500 Bohr. In FIG. H(b), the length L of the SNSNS system is 
8000 Bohr and the lengths of the N parts are again kept at 1000 Bohr. This makes the 
middle superconductor 6000 Bohr in length, which is longer than the BCS coherence length. 
In both figures the gap of the middle superconductor is 0.25A. Thus, there is effectively 
complete decoupling in FIG ^(h). In addition, one sees that the (2, 2) peaks are lying closer 
to one another for the longer system, which is in line with the behavior of the energies. 

IV. SELF-CONSISTENT CALCULATION OF THE GAP 

In the preceding discussions, the gap function A used in the calculation is steplike, that 
is, it has a finite constant value in the superconducting part and is zero in the normal part. 
At the interfaces, it has a step discontinuity. This configuration is shown in FIG. ^ So in 
the calculation of the bound states in Sec. |IT1|, the proximity effect is not taken into account. 
However we want to demonstrate in this section that our formalism makes it possible to show 
that the actual gap function is not steplike, but decreases smoothly towards the interface and 
abruptly goes down to zero in the normal metallic layer. This proximity effect is studied in 
the present section. In addition, the actual gap function will be calculated self-consistently. 

The self-consistency condition is given by 

A(x) = -yF(r,r,0+) (27) 

in which F{r, r, 0"*") is the well-known anomalous Green's function given by the upper right 
element of the original (r, r)-dependent matrix Green's function. Taking into account the 
expansion of the matrix Green's function over the Matsubara frequencies and the transverse 
wave vectors, we obtain 
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A(a:) = - F{x,x,ky,kz,iujn) (28) 

y ^ UJfi ■i^y i^z 

where F{x,x,ky,kz,iuJn) is the upper right element of the matrix Green's function in Eq. 
(pop. The summation over the Matsubara frequencies is restricted by the Debye temperature 
60 according to the formula, 

ksOD = ^nma. = nynaxT^ksT. (29) 

The summation over the transverse wave vector is over all positive values of ky and kz 
according to Eq. @). Evaluation of the summation takes much computer time, so we resort 
to some approximations. Since the cross section is a square, that excellent 
approximation to reduce the number of terms is to partition the transverse fc^-plane by 



concentric circles with radius fc_L = ^/c^ + k\ and a bk^ = j^. The number of allowed values 
of ky and k^ in each ring is given by its surface divided by the density of the transverse 
k states (^-^^ . Another approximation can be made by noting that the terms involving 
k± » ^JJi do not have significant contributions to the sum. This allows us to evaluate a 
finite number of terms instead of evaluating an infinite series. In our calculation we take 
fc±„aa: = 3^/1. In determining the coupling constant V according to the BCS-relation 

= 1.13 ujb , with A^(/i) = ^ (30) 

we use Tc = 1.2 K and ujo = 375 K for aluminum. We find V = 9.516 Ry. 

A. The Bar-Shaped Superconductor 

The matrix Green's function for a bar-shaped superconductor G^g{ X, X, ky^ k^^ '^^n) is given 
by Eq. (|n|). By substituting the upper right element of this Green's function in Eq. ( P5[ ) 
one straightforwardly obtains the following selfconsistency condition for the homogeneous 
bar-shaped superconductor, 

^^""^ ^ ^S,'^. {'^^TWs ^ TC^s) ' ^^^^ 
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In calculating the gap self-consistently, we first assume an initial value of the gap. By 
substituting it on the right-hand side of Eq. ( |3TD we obtain a new value of the gap. Then 
this new value is again substituted on the right-hand side to get another new value. This 
procedure can be repeated until the difference between successive iterations is negligibly 
small. As shown in FIG. ^ the difference between the first ten iterations is still significant. 
After about 80 iterations, the value of the gap stabilizes to 1.106 x 10~^ Ry. The transverse 
length is 1000 Bohr and the temperature is 0.6 K. The initial value of the gap is 2.0 x 10~^ 
Ry. 



FIG. |10] shows the plot of the gap against the transverse length Lt for different tem- 
peratures. The number of iterations is 100. It can be seen that there are oscillations of 
the gap. The amplitudes of the oscillations decrease as the transverse dimension increases. 
These oscillations can be attributed to the discreteness of the transverse wave vector. As 
the transverse width increases, the transverse wave vector approaches the continuous regime 
which can be gauged from the gap becoming closer to its bulk value obtained by integrating 
instead of summing over the transverse wave vectors. In the figure, we show the bulk value 
at 0.2 K. Another interesting thing which can be seen in the figure is the suppression of 
superconductivity for narrower transverse dimensions. We notice that as the temperature 
increases the onset of the suppression of superconductivity occurs at higher values of Lt. 



In FIG. |Tl] we show the temperature variation of the gap for different transverse widths 
Lt- A residual value of the gap for the bulk superconductor can still be observed beyond 
the critical temperature (1.2 K for aluminum). We also notice small- amplitude oscillations 
of the gap at higher temperatures, which only show up in the curves for larger transverse 
dimensions. For smaller transverse dimensions, the oscillations are suppressed. These oscil- 
lations are due to the cut-off in the summation over the Matsubara frequencies. For lower 
temperatures, the cut-off value n^nax in Eq- ( |29D becomes very large and the results are no 
longer sensitive to it. 
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B. The NS System 



The matrix Green's function appropriate for a normal metal-superconductor system is 
given in Eqs. (|l^) and (|l^. The first term on the right-hand side of Eq. ([T6| ) is the matrix 
Green's function for the bar superconductor. The second term contains the elements of the 
t-matrix, which take into account the scattering of the quasiparticles at the interface. The 
latter acts as a perturbing term to the former and is therefore responsible for the spatial 
variations of the gap and the pair amplitude at the vicinity of the interface. To calculate 
the spatial variation of the gap according to Eq. we substitute the upper right element 
of the matrix Green's function, Eq. ([16|), using the value of the self-consistent gap for the 
bar-shaped superconductor. 

Before presenting our results, we want to mention a computational problem which has to 
be solved. Due to the explicit presence of the position xj in Eq. ([T8|) , the matrices involved 
have alternating columns of very small and very large values which the computer cannot 
handle anymore. However this problem is not intrinsic to the formalism and can be solved 
by appropriate rescaling. By defining the matrix i 

^ ^'^-K.-^i^^^i'^'^'K:.^^ (32) 

the Green's function (|TBp, combined with Eq. (0), obtains the form 

a 

+ E - ^.oci^'^'C'Ti^' - ^j) (33) 

era' 

in which only position differences occur. The rescaled matrix i is determined by the equation 

y: z/c,<;(o)c°:7' = -^'<:;'(o), (34) 



found straightforwardly from the original equation 

FIG. p!2|(a) shows the spatial variation of the gap near the interface of a normal metal- 
superconductor system at T = 0.6 K for different transverse widths. At a distance of 30000 
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Bohr from the interface, which is about six times the coherence length 4500 Bohr), the gap 
is shghtly smaller than the one obtained for the bar-shaped superconductor. The differences 
are about 7.79%, 6.57%, and 2.93% of their bar value for Lt = 99.8514, Lt = 100, and 1000 
Bohr, respectively. At 4000 Bohr, which is of the order of the coherence length from the 
interface, the differences are about 21.7%, 20.44%, and 14.71%, respectively. These figures 
lead to the inevitable conclusion that for larger transverse dimensions, the proximity effect is 
less pronounced than for smaller ones. This may be due to the fact that for small transverse 
dimensions the superconductivity tends to be suppressed. We find that the difference seen 
for the two smaller widths can be attributed to the oscillations in the self-consistent gap 



shown in FIG. [1^. We have not seen a special influence of the fact that the smaller width 
is a critical one. FIG. ^{h) shows the corresponding pair amplitude defined by (see 
Eq. (pT])) which has a finite value in the normal metallic region near the NS interface but it 
decays in the inner region of the normal metal. 

C. The SNS System 

In the superconductor-normal metal-superconductor system, there are two interfaces 
which we designate as Xi and X2- The t-matrices must be determined at these interfaces so 
that we can evaluate the T-matrix. The unpleasant singularities due to the explicit presence 
of the position of the interfaces in Eq. ([18|) can be removed by using Eq. (^21) resulting to 
a corresponding transformation of the T-matrix given by 

f;/;;;r' = e*--'^^i"^T;/;'<r'e*"'"''^^v"^'. (35) 

By implementing this, the multiple-scattering Green's function ( ]20|) and Eq. (pT]), which 
determines the T-matrix in terms of the t-matrices, can be modified straightforwardly. 

The steps in calculating the gap and the pair amplitude self-consistently are the same 
as in Sec. [IV B|. FIG. IT3| shows the self-consistent gap function and the pair amplitude for 



different transverse widths. The center of the system is at x = 0. The spatial variation of 
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the gap near the interfaces is clearly shown. The proximity effect is stronger than for the 



NS system shown in FIG. O. Whereas in FIG. 12a for Lf = 1000 Bohr the gap at a distance 



of 6000 Bohr from the interface has a value of 9.753 x 10^ Ry, in FIG. |T3|a it has already 
increased to 1.074 x 10~^ Ry, which is much closer to the bulk value of 1.106 x 10~^ Ry. In 
FIG. |I3]b the pair amplitude in the N region does not decrease below a value of 5 x 10^^. In 



FIG. [12| b it decreases to zero and the value at 2000 Bohr from the NS interface has already 
decreased to 2.784 x 10"''. In FIG. |14| we show the gap functions for different values of L. 



It can be seen that for larger L the gap function is lesser in magnitude. This is another 
manifestation of the proximity effect. 

D. The SNSNS System 

The extension of the steps outlined in Sec. [IV C| leading to the T-matrix and the matrix 
Green's function for a SNSNS system is straightforward. In this case we have to consider 
four interfaces but the procedure is basically the same. In the outer superconductors, we 
again use the self-consistent value of the gap for the superconducting bar we have calculated 



in Sec. [IV A| . In calculating the gap in the middle superconductor self-consistently, we apply 
the recipe introduced by Tanaka and Tsukada |[l[] . We start the iteration procedure by taking 
half of the gap in the outer superconductors for the gap in the middle. By substituting these 
values in the righthand side of Eq. (p8[) , the middle a;-dependent gap function is obtained 
and its spatial average is determined. This average value is used as the new value of the 
gap in the middle superconductor. The process is repeated until the difference of these 
average values between successive iterations is negligibly small. To determine the gap's 
spatial variation, the self-consistent average value of the gap in the middle superconductor 
and the bar value in the outer ones are substituted in the right-hand side of Eq. 



FIG. shows the gap and the pair amplitude of two SNSNS systems, one with a length 
Ln = 3000 Bohr for the normal metallic parts and the other one with a larger length 
-Ltv = 5000 Bohr. First note that the gap of the middle superconductor is smaller than the 
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gap of the outer superconductors although the metalhc parameters for these superconductors 
are taken equal. This supports our choice of a /i value smaller than 1 in Sec. |T|, see FIG. |^. 
Further note that the gap and the pair amplitude are higher for the system with narrower 
normal-metal part, which is just another manifestation of the proximity effect. Although 
by definition the gap is zero in the normal-metal parts, the pair amplitude is different from 
zero. This is also a manifestation of the proximity effect. If the width of the normal-metal 
layers had been much larger than the sum of the coherence lengths of the outer and middle 
superconductors, the pair amplitude would have been zero there except in those regions near 
the interfaces. 

V. CONCLUSIONS AND FUTURE PROSPECTS 

The multiple-scattering Green's function formalism developed by Koperdraad has 
been apphed to determine the Andreev bound states in SNS and SNSNS junctions through 
the local density of states and to calculate the superconducting gap function self-consistently. 
We have shown that for transverse junction widths tuned such that the motion in the longi- 
tudinal direction is very slow, the Andreev approximation breaks down. For these transverse 
dimensions, the highest mode is supported by many bound states whose degeneracy is lifted 
when no approximation is applied. The thickness used for the normal metallic layers is 
chosen such that the lower modes are each supported by only one nondegenerate state. Re- 
sults for the self-consistent gap functions exhibit various features of the proximity effect. 
Furthermore, our results show that for small transverse dimensions, superconductivity is 
suppressed. 

The formalism is applicable to systems of an arbitrary number of layers. In addition, 
it allows for the calculation of the supercurrents through such junctions, to which self- 
consistent gap functions are necessary [0,0. This will be investigated in the near future. 
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CAPTION FOR FIGURES 



FIG. ^ The geometry of the bar-shaped superconductor under consideration. It has 
finite transverse dimension and is infinite longitudinally. A rectangular cross- section is 
shown whose dimensions are Ly and L^. 

FIG. ^ A two-dimensional picture of the systems studied. In the upper panel the vertical 
direction stands for a transverse direction. In the lower panels the potential is shown, which 
is zero in the normal part(s), and proportional to A in the superconducting parts. 

FIG. H Plot of the LDOS against E/A for a SNS system in which Lt = 13 Bohr. 

FIG. § Plot of the LDOS against E/A for x = 1000 Bohr, Lt = 12.5676 Bohr, L=4000 
Bohr, and A = 0.0001 Ry. 

FIG. |. Plot of the LDOS for a SNS system against E/Aatx = 1000 Bohr (solid curve) 
and X = 1500 Bohr (dashed curve). The transverse dimension is Lt = 12.5676 Bohr and the 
length of the normal metal is L = 4000 Bohr. 

FIG. |. Plot of the LDOS against E/A for a SNS system in which Lt = 100 Bohr (solid 
curve) and Lt = 99.8514 Bohr (solid curve). The length of the normal-metal part is 4000 
Bohr. 

FIG. 0. Plot of the LDOS at a; = -1500 Bohr against E/A for a SNSNS system in which 
Lt = 12.5676 Bohr for (a) h = 0.25 (b) h = 0.5 (c) h = 0.75 and (d) h = l. 

FIG. |. The LDOS against E/A for a SNSNS system in which Lt = 12.5676 Bohr. 
In (a) L = 6000 Bohr, the length of the middle superconductor is 4000 Bohr and the 
LDOS is calculated at x = -2500 Bohr. In (b) L = 8000 Bohr, the length of the middle 
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superconductor is 6000 Bohr and the LDOS is calculated at a; = —3500 Bohr. The gap of 
the middle superconductor is 0.25A. All unlabeled LDOS peaks belong to the mode (2,2). 
The calculation is done in the Andreev approximation. 



FIG. ^. The gap against the number of iterations for a bar-shaped superconductor. Note 
that the value of the gap stabilizes as the number of iterations increases. For the system 
considered, Lt = 1000 Bohr and T = 0.6 K. 

FIG. |10|. Plot of the self-consistent gap function against the transverse length Lt for a 
bar-shaped superconductor at different temperatures. The number of iterations is 100. 



FIG. |TT]. The temperature variation of the self-consistent gap for different transverse 
widths. 

FIG. [l^. (a) The gap and (b) the pair amplitude against the distance from the interface 
of a NS system at T = 0.6 K. The interface is chosen at x = 0. 

FIG. |13|. (a) The gap function and (b) the pair amplitude of a SNS system against 
the distance from the middle of the system chosen at x = 0. The interfaces are located at 
X = ±2000. 



FIG. 0. The gap function of a SNS system against the distance from the middle of the 
system chosen at x = for different values of L. The interfaces are located at ±2000 for 
L = 4000 Bohr and at ±4000 for L = 8000 Bohr. 

FIG. |15[ (a) The gap function and (b) the pair amplitude against the distance from 
the middle of a SNSNS system, which is chosen at x = 0. The interfaces are chosen at 
X = ±2000, ±5000 for Ln = 3000 Bohr (sohd curve) and x = ±2000, ±7000 for = 5000 
Bohr (dotted curve). The transverse width is Lj = 100 Bohr. 
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FIGURES 



FIG. 1. The geometry of the bar-shaped superconductor under consideration. It has finite 
transverse dimension and is infinite longitudinally. A rectangular cross-section is shown whose 
dimensions are Ly and L^. 
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FIG. 2. A two-dimensional picture of the systems studied. In the upper panel the vertical 
direction stands for a transverse direction. In the lower panels the potential is shown, which is zero 
in the normal part(s), and proportional to A in the superconducting parts. 
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FIG. 3. Plot of the LDOS against Ejt^ for a SNS system in which Lt = 13 Bohr. 
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FIG. 4. Plot of the LDOS against E/A for x = 1000 Bohr, Lt = 12.5676 Bohr, L=4000 Bohr, 
and A = 0.0001 Ry. 
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FIG. 6. Plot of the LDOS against £;/A for a SNS system in which Lt = 100 Bohr (dotted 
curve) and Lt = 99.8514 Bohr (sohd curve). The length of the normal-metal part is 4000 Bohr. 
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FIG. 7. Plot of the LDOS at x = -1500 Bohr against E/A for a SNSNS system in which 
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FIG. 8. The LDOS against ^/A for a SNSNS system in which Lt = 12.5676 Bohr. In (a) 
L = 6000 Bohr, the length of the middle superconductor is 4000 Bohr and the LDOS is calculated 
at X = —2500 Bohr. In (b) L = 8000 Bohr, the length of the middle superconductor is 6000 Bohr 
and the LDOS is calculated at x = —3500 Bohr. The gap of the middle superconductor is 0.25A. 
All unlabeled LDOS peaks belong to the mode (2,2). The calculation is done in the Andreev 
approximation. 
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FIG. 9. The gap against the number of iterations for a bar-shaped superconductor. Note that 
the value of the gap stabihzes as the number of iterations increases. For the system considered, 
Lt = 1000 Bohr and T = 0.6 K. 
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FIG. 10. Plot of the self-consistent gap function against the transverse length Lt for a 
bar-shaped superconductor at different temperatures. The number of iterations is 100. 
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FIG. 11. The temperature variation of the self-consistent gap for different transverse widths. 



29 




1.2x10* 




X (distance from the middle of tfie system in Bohr) 



(b) 



FIG. 12. (a) The gap and (b) the pair amplitude against the distance from the interface of a 
NS system at T = 0.6 K. The interface is chosen at a; = 0. 
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FIG. 13. (a) The gap function and (b) the pair ampHtude of a SNS system against the distance 
from the middle of the system chosen at a; = 0. The interfaces are located at x = ±2000. 
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FIG. 14. The gap function of a SNS system against the distance from the middle of the system 
chosen at x = for different values of L. The interfaces are located at ±2000 for L = 4000 Bohr 
and at ±4000 for L = 8000 Bohr. 
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FIG. 15. (a) The gap function and (b) the pair amphtude against the distance from the middle 
of a SNSNS system, which is chosen at a; = 0. The interfaces are chosen at x = ±2000, ±5000 for 
Ln = 3000 Bohr (sohd curve) and x = ±2000, ±7000 for L;v = 5000 Bohr (dotted curve). The 
transverse "width is = 100 Bohr. 
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